A remedy for constraint growth in Numerical Relativity 
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Rapid growth of constraints is often observed in free evolutions of highly gravitating systems. 
To alleviate this problem we investigate the effect of adding spatial derivatives of the constraints 
to the right hand side of the evolution equations, and we look at how this affects the character 
of the system and the treatment of boundaries. We apply this technique to two formulations of 
Maxwell's equations, the so-called fat Maxwell and the Knapp-Walker-Baumgarte systems, and ob- 
tain mixed hyperbolic-parabolic problems in which high frequency constraint violations are damped. 
Constraint-preserving boundary conditions amount to imposing Dirichlet boundary conditions on 
constraint variables, which translate into Neumann-like boundary conditions for the main variables. 
The success of the numerical tests presented in this work suggests that this remedy may bring 
benefits to fully nonlinear simulations of General Relativity. 



I. INTRODUCTION 

The 3 + 1 decomposition of Einstein equations gives 
rise to a system of evolution equations and a set of con- 
straints. Although it is true that, at the analytical level 
and in the absence of boundaries, if the constraints are 
satisfied on an initial spacelike hypcrsurface, they remain 
satisfied at later times by virtue of the evolution equa- 
tions, this is not the case in numerical simulations. If 
the continuum initial data satisfy the constraints exactly, 
then they will only satisfy the discretized constraints up 
to truncation error. In many cases, initial errors grow ex- 
ponentially, and sometimes even more rapidly. Whether 
the violent growth of the constraints is causing the crash 
of numerical codes, or is simply the effect of a more gen- 
eral growth in the error, is not altogether clear. However, 
the hope is that, if the constraint growth can be pre- 
vented or at least alleviated, the runtime of numerical 
simulations may increase. 

In recent years a number of techniques have been pro- 
posed aimed at damping the constraints or, at least, con- 
trolling their growth. These are the A-system the fine 
tuning of parameters multiplying the constraints that are 
added to the evolution equations Q, the dynamical ad- 
justment of such parameters and the addition of vari- 
ational derivatives of a constraint energy to the evolution 
equations The method that we wish to investigate 
has remarkable similarities with the latter. When time- 
like boundaries are present all of these methods require 
constraint-preserving boundary conditions to prevent the 
injection of constraint violations. 

We begin with a simple observation. Consider the first 
order scalar equation on the real line, 



Ui 



(1) 



where : 



-1, supplemented with smooth initial data. If 



we make the ansatz u(t, x) = e %UJX u(t, u>), to e K, and sub- 
stitute it into P| ; we obtain u t = —OJU. The fact that the 



solution u(t, x) = e 1 ' 



u(0,u>), allows for unbounded 



in a stable way. Now we modify our problem and assume 
that the constraint 



C = u x + iu = 



(2) 



has to hold. Clearly, if u satisfies Eq. JJJ, then Ct = 
iC x , i.e., the constraint propagation problem is also ill- 
posed. A crucial point is that we are allowed to add 
constraints and/or derivatives of constraints to the right 
hand side of the evolution equations without affecting the 
space of physical solutions (i.e., solutions that satisfy the 
constraints). For example, adding a spatial derivative of 
the constraint to the right hand side of Eq. JTJ leads to 



u t = iu x + C x = 2iu x + u xx 



(3) 



growth |2 
posed and 



shows that the initial value problem is ill— 
therefore, cannot be consistently discretized 



This equation is parabolic and, most importantly, gives 
rise to a well-posed initial value problem for both the 
main and the constraint propagation systems. Not 
only does the problem no longer suffer from unbounded 
growth of the constraints, but the modification also 
damps the high frequency constraint violations. Fur- 
thermore, in the presence of boundaries, we can impose 
homogeneous Dirichlet boundary conditions on the con- 
straint, C = 0, which, when translated in terms of the 
main variables become of Neumann type, u x = —iu. This 
example, although trivial, clearly illustrates the poten- 
tial benefits of dissipation. It allowed us to convert an 
ill-posed initial value problem into a well-posed initial- 
boundary value problem |21| . 

In this paper we investigate the effect of adding spa- 
tial derivatives of the constraints to the evolution equa- 
tions for two formulations of Maxwell's equations, the 
so-called fat Maxwell [j| (first order) and the Knapp- 
Walker-Baumgarte (KWB) 0| (first order in time, sec- 
ond in space) systems. As expected, the modification 
destroys the hyperbolic character of the systems. How- 
ever, if done carefully, it can lead to mixed hyperbolic- 
parabolic systems both for the evolution equations and 
the evolution of the constraints, such that the high fre- 
quency components of the constraints are damped. 

In the presence of timclikc boundaries the systems that 
we analyze are such that constraint-preserving boundary 
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conditions amount to imposing homogeneous Dirichlet 
boundary data on the constraints, which translate into 
Neumann-like boundary conditions for the main system. 
Compared to standard constraint-preserving boundary 
conditions for hyperbolic problems yj S l1 OH [HI > those 
that we obtain can be much simpler and possibly more 
robust. For example, if the constraint propagation sys- 
tem is strongly parabolic, then at the boundary we can 
impose that all constraints vanish, which is less involved 
than computing the incoming characteristic constraints 
and, as usually done, trading the normal with time and 
tangential derivatives using the evolution equations. Fur- 
ther, the modified problems allow for a greater number of 
constraint components to be set to zero. Unfortunately, 
— and this is presumably the most regrettable feature of 
this technique — due to the loss of finite speed of prop- 
agation, the time step needs to be proportional to the 
square of the mesh spacing in order to achieve numerical 
stability (and convergence) with explicit finite difference 
schemes. 

Parabolic systems have received little attention in nu- 
merical relativity with some notable exceptions [l2^ . in 
which the aim was primarily to promote elliptic gauge 
conditions to parabolic equations. Although the tech- 
nique of adding spatial derivatives of the constraints has 
been previously employed by Yoneda and Shinkai fl3l ] 
and by Fiske |4|, the analysis of the character and the 
study of appropriate boundary conditions for the modi- 
fied systems, to our knowledge, is new. 

We recall some basic definitions that will be used 
throughout the paper. Consider the systems of partial 
differential equations in one space dimension 

d t u = Ad 2 x ii + Bd x u + Cu, (4) 
d t v = Dd x v + Ev, (5) 

where A, B, C, D, and E are constant square matrices. 
System (0J is said to be strongly parabolic if the eigenval- 
ues of A + A^ are positive. We call system (0 symmetric 
hyperbolic if D = D> . In the presence of boundaries and 
with appropriate initial and boundary conditions both 
systems admit a well-posed initial-boundary value prob- 
lem as does the coupled mixed hyperbolic-parabolic sys- 
tem [3 

„ / u \ ( Adl+Bd x + C Fd x + G\ ( u\ 
d \v) = \ Ld x + M Dd x + E){v) < 6 ) 

for any choice of the coupling matrices F, G, L, and 
M . Notice that the character of system © is entirely 
determined by the two matrices A and D. 

The paper is organized as follows. In Sec. [H] we in- 
troduce the fat Maxwell system, appropriately add spa- 
tial derivatives of constraints and analyze the resulting 
initial-boundary value problem. A similar analysis is 
then repeated for the KWB system in Sec. lIIII In Scc. lIVI 
after describing the discretization of the system and, par- 
ticularly, of the boundary conditions, we summarize the 
results of our numerical experiments. We conclude with 



some final remarks in Sec.0 In appendix lAl we discuss 
an energy conserving discretization for the wave equation 
in second order form, allowing for Sommerfeld boundary 
conditions. This discretization is used for the hyperbolic 
sector of the KWB system. 

We use units where c = /io = £o = 1 and adopt the 
Einstein summation convention. 



II. THE FAT MAXWELL SYSTEM 

The vacuum Maxwell equations on a Minkowski back- 
ground in Cartesian coordinates are given by 

d t Ei = <,,,i) ! n" . (7) 

dtBt = -e ijk d j E k , (8) 

d k E k = 0, (9) 

d k B k = 0. (10) 

In terms of the potentials, <p and Ai, the electric and 
magnetic fields are 

Et = -d t Ai-d i( j>, 

Bi = ei jk d j A k . 

If we introduce the variables Dij = diAj, we can rewrite 
Maxwell's equations in the form [f| 

d t Ei = d k (D ik -D ki ), 

d t Dij = -diEj - d % dj4> , 

C = d k E k = , 

Cijk = diD jk - djD ik = , 

where we used the identity e^e 11 ™ = <5j<5™ — S™5 l k . This 
system has 12 variables that depend on the space-time 
coordinates and is subject to 10 constraints. One can 
readily verify that the constraints propagate trivially, 

d t C = 0, 
dtCijk = , 

and that the evolution system is weakly hyperbolic, i.e., 
it cannot be diagonalized. However, by appropriately 
adding constraints to the right hand side of the evolution 
equations 

dtEi - d k (D ik -D ki )+C ikk 

= -d k D kl + d x D kk , 

d t Dij = -diEj - didj(j) + SijC 

= -d i E j -d i d j <j> + 5 ij d k E k , 

one can ensure that both the main evolution system and 
the evolution of the constraint variables are symmetric 
hyperbolic. We note that this modification corresponds 
to setting 71 = 2 and 72 = 1 in the two parameter family 
of formulations of Lindblom et al. [a. 
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Any first order, linear, homogeneous, constant coeffi- 
cient, symmetric hyperbolic system with no lower order 
terms admits, by definition, a conserved energy. An en- 
ergy, in this context, is an integral over the spatial do- 
main, f2, of a positive definite quadratic form of the main 
variables. The energies 

E = \ [ (EfiEi + ) d 3 x , (11) 

Ec = \ j n {° 2 + \ CijkCi *) d3x ' ( 12 ) 

for example, are conserved up to contributions due to the 
boundary (the conservation of E requires the Coulomb 
gauge, = 0, which will henceforth be assumed). An 
immediate consequence of this is that, if = R 3 or fl = 
T 3 , the energy associated with an initial violation of the 
constraints will be preserved during evolution. We now 
show that by adding spatial derivatives of the constraints 
to the evolution equations it is possible to improve on 
this. 

The general idea is to add constraints and/or deriva- 
tives of constraints to the right hand side of the evolu- 
tion equations that generate decaying or damping terms 
in the constraint propagation system. Consider the fol- 
lowing modification of the main evolution system 

dtEi = ... + fidiC (13) 

= -d k D kl + d t D kk + »did k E k , 

dtDij = ... + f,d k C kij (14) 
= -diEj + S l3 d k E k + txd k {d k Dij - diD kj ) . 

The evolution of the constraint variables becomes 

d t C = ffCok + nffdiC, (15) 
d t C ijk = 2d [i S j]k C+2fxd li d l C mk . (16) 

By following this procedure, we have lost symmetric 
hyperbolicity. It is essential to check whether the mod- 
ified systems still admit a well-posed initial value prob- 
lem. We start by showing that both systems, (|13[l 114|) 
and p5 jl -l|16 p . satisfy a necessary condition for well- 
posedness, the Petrovskii condition |14|. We will show 
that for positive values of the parameter /i most Fourier 
components of the constraint variables are damped. We 
then give an energy estimate for the evolution of the con- 
straints systems and, as a byproduct, obtain constraint- 
preserving boundary conditions. Our attempts at de- 
riving an estimate for the main system have not been 
successful when boundary conditions consistent with the 
constraints are imposed. Nevertheless, we found the one 
dimensional reduction of the problem to be well-posed 
and therefore used it as a tool for obtaining appropriate 
boundary conditions for the main system. The rationale 
being that any boundary condition that makes the re- 
duced problem ill-posed would be unsuitable in higher 
dimensions. 



Let Ut = P(d x )u be a system of linear constant coef- 
ficient partial differential equations in d dimensions. As- 
sume for the moment that there are no boundaries. Ac- 
cording to the Petrovskii condition, it must be possible to 
bound the real parts of the eigenvalues of the symbol, or 
Fourier transform, of the differential operator, P(ico), by 
a real constant a, which is independent of u>. For scalar 
equations this condition is also sufficient. 

To compute the symbols of l(IS ty -{T2 |l and JTSJ-JT^I we 
use the substitution d k — > iuj k = i\w\u) k , where <l> k <l> k = 1. 
From the evolution equations we get 

d t E u = i\uj\D AA - huj 2 E u , 

d t E A = -i\u\D u A , 

d t D uu = 0, 

dtD uA = -i\u\E A , 

d t D Au> = -[lu?D Au} , 

d t D AA = 2i\oj\E u - imjj 2 D aa , 

d t D A ~B = -^ d Ab E ^ ' 

where u u = u k u k , u A = A k u k , A k B k = S AB , uj k A k 
0. and U£ B = u A b — \5 A bucc- The evolution of the 
constraints yields 

d t C = i\u>\C uAA - i^uj 2 C , 
d t C uA A = 2i\w\C - (jio 2 C uAA , 
d tC uAB = -fMUJ C uAb , 
dtC uA w = —V^CuAm j 
d t C AB k = 0. 

The eigenvalues of the Fourier transformed systems form 
a subset of 

{±i|w|, —fiuj 2 ± iV2\u>\, —^uj 2 , 0} . 

Hence the Petrovskii condition is satisfied provided that 
[i > 0. In particular, this shows that for [i > the high 
frequency components of most constraints arc rapidly 
damped. 

We now calculate the time derivative of (|12fl . This is 
given by 

-E c = / {CC nAA + vCdnC + fxC nAk diC lAk ) d 2 S 
dt Jan 

-M / (d'CdiC + diC^^dj^^x. (17) 
Jn 

Clearly, the volume integral gives a non-positive con- 
tribution to the time derivative of the energy |22l| and, 
therefore, by controlling the boundary term one can en- 
sure that the energy associated with the constraints is 
non-increasing. Boundary conditions that can be used 
for this purpose are 

C = 0, O nAk = 0. (18) 

When these 7 conditions are enforced we obtain the in- 
equality 

^-E c = -n f (d l Cd z C + diC ijk d l C ljk ) d 3 x<0. 
dt Jn 
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We now turn to the main system. As mentioned earlier, 
in the presence of boundaries we were not able to obtain 
an energy estimate using conditions p8|l. We therefore 
look at simple necessary conditions for well-posedness for 
the quarter space problem. We assume that the bound- 
ary is located at x = and that £1 = {(x, y, z) £ M. 3 \x > 
0}. If the variables depend only on x and t, the system 
becomes 

d t Ei = +d x D AA + l xd 2 x E 1 , 
d t E A = -d x Di A , 
d t Dn = 0, 
d t D 1A = -d x E A , 
d t D A i = [id 2 x D A1 , 
d t D AB = 5 AB d x E 1 + [id x D A B , 

where A = 2, 3 and B = 2, 3. No matter what 
boundary conditions one uses for the main system, they 
should be such that the one dimensional reduction of 
the problem is well-posed. The equations above form 
a mixed hyperbolic-parabolic system. This can be seen 
by identifying u = {Ex,D Ak ,} and v = {E A ,Di k }, 
and comparing with Eq. JBJ. The 7 variables belong- 
ing to the parabolic sector require a boundary condition 
each. Eq. I|18|) provides us with the correct number of 
Neumann-like boundary conditions. We are free to spec- 
ify maximally dissipativc boundary conditions to the 5 
variables of the hyperbolic sector, 

{w Ain - S AB W B out)\x=0 = 9A , (19) 

where S is a sufficiently small 2x2 matrix, g A is freely 
specifiable boundary data and 



where the fields A;, Ei and T are subject to the con- 
straints 



1 



W A ir 



W Aout 



= -j={E A + D 1A ), 



T2 iEA 



Dxa), 



are the ingoing and outgoing characteristic variables. 
This gives a total of 9 boundary conditions for the main 
evolution system, 7 of which ensure constraint preserva- 
tion m. 

A description of the discretization of the initial- 
boundary value problem and the results of the numerical 
experiments are given in Sec. II VI 



III. THE KWB SYSTEM 

In Ref. Knapp, Walker, and Baumgarte, introduced 
a formulation of Maxwell's equations which shares some 
features with the Baumgarte-Shapiro-Shibata-Nakamura 
formulation of Einstein equations |l5lll6j . In vacuum and 
on a Minkowski background in Cartesian coordinates the 
equations take the form 



d t Ai = -Ei - di<j) , 

dtBk = n,n-A, ■ o, i . 

d t r = ~djd j (t>. 



(20) 
(21) 
(22) 



C E = d l E t = , 
CV = d;A l - T = . 



(23) 
(24) 



The system contains second spatial derivatives and is 
symmetric hyperbolic 01 ■ As the authors of [|| pointed 
out, in this system the constraints propagate according 
to the wave equation 



d t C E = -8^ Or 
dtCr = —Ce , 



(25) 
(26) 



which is to be contrasted with the static constraint evo- 
lution of the original Maxwell system, Eqs. (fzfl — (| 1 C)|) . 

The constraint system (|25|) - (|26|l conserves the follow- 
ing energy 

Ec = \ f(C 2 E + a i C r ^C r )d 3 x . (27) 

If one sets <j) — and modifies the evolution equations 
according to 

d t At = ... + Lid t C r (28) 

= -E i + fxd i d j A j -fj,diT, 
d t Et = ... + fidiC E (29) 

= ,),;>* A, + diT + ,!<),<)■ e, , 

d t T = XCr = XdiA 1 - Ar . (30) 

which can be identified with Eqs. (22)-(23) of Fiske 0], 
the evolution of the constraints becomes 

d t C E = -d^iCr + n&BiCE , (31) 
d t C v = -CE + ^d l d l C v -XC v . (32) 

The evolution equations were modified so that 1131(1 and 
would acquire damping and friction terms. 
Rather than estimating the energy (|27|l we observe 
that, if /x > 1/2, the constraint propagation system is 
strongly parabolic and an energy estimate for the con- 
straint system can be obtained directly in L^. The time 
derivative of 



Ec = \ f (C 2 E + C 2 ) d 3 



(33) 



d_ 
~dt 



-E c = j (-C E d n Cr + vCEd n CE + Crd n Cr)d 2 S 
Jan 

+ [ (WCediCr - fid l C E diC E - 0C T diC r ) d 3 x 
Jn 

+ ! {-C T C E - XCl) d 3 x . (34) 

JQ. 

The use of Dirichlet boundary conditions, 

C E = C T = , (35) 
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ensures that the boundary term of (|34[) vanishes. Fur- 
thermore, the volume term containing spatial derivatives 
is non-positive. Thus, an energy estimate can be ob- 
tained |24| . 

We now place a boundary at x — so that SI = 
{(x,y,z) £ M. 3 \x > 0}. By ignoring any y and z de- 
pendence the hyperbolic-parabolic character of the sys- 
tem becomes evident (the introduction of the auxiliary 
variable Wb = 9 x Ab may be helpful for this analysis) 



d t A 1 = 


-E 1 + [idlA x - 


- ^d x T , 


d t A B = 


-E B , 




d t E x = 


-dlA x + d x T - 


- ildlEy 


d t E B = 


-d*A B , 




d t v = 


\d x Ai - XT . 





For /i > 1/2 we can identify a parabolic sector {A\, E{\ 
and a hyperbolic sector {Ab 1 Eb, T}. Dirichlet boundary 
conditions on the constraints, Ce = 0, Cr = 0, corre- 
spond to Neumann-like boundary conditions on the main 
variables of the parabolic sector. Maximally dissipativc 
boundary conditions can be imposed on the variables of 
the hyperbolic sector, 

(d x A B + E B )\ x =o = . (36) 

Numerical experiments with the modified KWB sys- 
tem indicate that the discrete approximation of the 
initial-boundary value problem is stable. More details 
are given in the next Section. In Appendix^wc analyze 
the accuracy and stability of discrctized maximally dissi- 
pative boundary conditions, such as Eq. (|36|) . for the one 
dimensional wave equation in second order form. 

IV. NUMERICAL IMPLEMENTATION 

To discretize our problems we use the method of lines, 
whereby space is discretized while time is left continuous. 
We then integrate the resulting system of ordinary differ- 
ential equations using 4th order Runge-Kutta. We take 
our domain to be f2 = [0, l] 3 and use periodic bound- 
ary conditions in y and assume no dependence on z. We 
introduce a grid such that the points (0,j, k) belong to 
the boundary, x = 0. In the interior we approximate the 
first spatial derivatives with centered differencing opera- 
tors and the second spatial derivatives according to 

d id mu — > <f D +^ lD -x lU ijk I = m 

X ' ^ \ II,, D 0x ' ' "' 

where hD + Uj = itj+i — Uj, hD^uj = Uj — Uj—i, and 
2Dq = D + + D-. Thus, the discretization in the interior 
is second order accurate. 

The discretization of the boundary conditions is based 
on algorithms for which convergence proofs are available, 
at least for simple toy model problems 01- For example, 
one can show that for the one-way advective equation, 



Ut = u x , x > 0, one can use one-sided differencing at the 
boundary, = D + uq. For the heat equation, u t = u xx , 
x > 0, with Neumann boundary condition u x (t,0) = 
0, one can use the approximation DoUa(t) = which, 
when combined with the evolution equation at i = 
becomes = j-D+uq- Both cases give global second 
order convergence. 

To simplify the numerical implementation of the 
boundary conditions we found it convenient to introduce 
a ghost zone consisting of an extra grid point, i = — 1, 
to the left of the boundary, so that derivatives could be 
computed at i = as if it were an interior point. In the 
fat Maxwell case we use the following prescription to pop- 
ulate the ghost zone: we use the Neumann-like boundary 
conditions l|18|l and use second order extrapolation for 
the remaining fields, {Ea, Dik}. Explicitly, 

= Do x Ei-ojk + Do y E2- t ojk + Do z E3- t ojk , 

= Do x DAm;Ojk — -DoA-Dlm;Ojfe , 

= D +x D- x EA : ojk , 
= D+ x D_ x Di m -ojk ■ 

Note that when combined with a second order centered 
difference at i = 0, second order extrapolation is equiv- 
alent to first order one-sided differencing. Eq. (|19|l with 
Sab = 9 A = is then used to prescribe the incoming 
fields of the hyperbolic part. The new variables are com- 
puted by solving the system 

£ (uew) +i) (new) = Q 

p(new) n (new) _ ^(old) n (old) 
A ~ U 1A A ~ U \A 

In the KWB case the boundary treatment is slightly 
different due to the hyperbolic role that some second 
spatial derivatives play in the evolution equations. The 
theory of difference approximations for second order hy- 
perbolic systems is much less developed than the cor- 
respondent one for first order systems (for some re- 
cent developments see Ref. 0). Based on the energy 
method for semi-discrete systems, in Appendix we 
give a proof of second order convergence for a particu- 
lar discretization of the one dimensional wave equation, 
Utt = Uxx, x > 0, with Sommcrfcld boundary conditions, 
Ut(t,0) = u x (t,0). This result is used in the discretiza- 
tion of boundary condition l|36|l . 

As in the fat Maxwell case, we introduce a ghost zone 
and populate it using the following equations 

= Do x Ai-ojk + Do y A2-fijk + Do z A3-Qjk — Tojk , 

= DoxAsfljk + Es-0,jk , 

= D 0x Ei- jk + DoyE2-,ojk + D 0z E 3;0 jk , 
= D +x D- x EB;0jk , 
= D+ x D- x Tojk ■ 

The first and third conditions correspond to the vanish- 
ing of the constraints, Eq. I|35|) . The second condition 
is the discretization of Eq. I|3t)|l . done according to the 
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FIG. 1: To test the stability and accuracy of our discretization 
we perform a grid refinement test. The domain is Q = [0, l] 2 
with no dependence in the z direction and periodic boundary 
conditions in the y direction. At x — and i = lwe use the 
discrete boundary conditions described in Sec. llVl For both 
the fat Maxwell and KWB systems we plot Q = log 2 (||uh — 
w e ||/||it/i/2 — u e\\) where Uh and Uh/2 represent the numerical 
solutions computed with a mesh spacing of h = 1/100 and 
Me the exact solution, which is given by a wave travelling in 
an oblique direction with respect to the boundary. For the 
fat Maxwell system it is obtained from the potentials <f> — 
and A = i sin[7r(a; + 2y + VSt)](0, 0, 1), while for the KWB 

system from = and A = i sin[7r(a; + 2y + \/Zt)](2, -1, 0). 
In the fat Maxwell case the time step, k, is chosen so that 
A/ik/h 2 = 1, where h is the spatial mesh spacing, and we set 
H = 1/10. In the KWB case we used 2/ik/h 2 = 1, fx = 1 
and A = 10. Values of Q close to 2 indicate second order 
convergence for the main variables. 



prescription (|A7|) . The last two conditions correspond to 
second order extrapolation. 

We treat the x = 1 boundary in a similar way. Having 
verified that our discretizations are second order conver- 
gent for both the fat Maxwell and the KWB systems, see 
Fig. ^ we test our schemes with constraint violating ini- 
tial data corresponding to a sufficiently smooth pulse of 
compact support to all fields and perform runs at three 
different resolutions. Given that in the fat Maxwell sys- 
tem there is no restriction on the value of fj,, we choose 
something reasonably small, 1/10, so that we can take 
larger time steps. In the KWB case, instead, we must 
choose fx > 1/2 otherwise it would not be clear that one 
could use Ce = Cr = as boundary conditions for the 
constraints. We opt for (j, = 1. The other parameter that 
is present in the KWB formulation is A. This is subject 
to no restriction, does not affect the choice of time step 
(for sufficiently high resolutions), and, if positive, con- 
tributes to the damping. We choose A = 10. Not only do 
we observe decay of the constraints, as shown in Fig. 
but we also see no signs of instability in the main system. 
Finally, in Fig. 01 we compare the constraint decay of the 
modified systems with that of the original ones. 




FIG. 2: We plot the value of E^ 2 , Eqs. fT^ and l|27> . for the 
fat Maxwell system (top) and the KWB system (bottom) for 
three different resolutions. The curves are labeled according 
to the number of grid points used in each spatial direction. We 
give constraint violating initial data consisting of a sufficiently 
smooth pulse of compact support to all fields. In both cases 
the domain is f2 = [0, l] 2 , having dropped the dependence in 
the z direction, and we impose periodic boundary conditions 
in the y direction. At x = and x = 1 we use the boundary 
conditions outlined in Sec. II VI The values of the parameters 
and the size of the time step k are chosen as in Fig. The 
noise in the curves of the KWB case are likely due to the fact 
that the energy being monitored, Eq. 1271 . contains spatial 
derivatives which are discretized using Do operators. 



V. CONCLUDING REMARKS 

In this work we investigated the effect of adding spa- 
tial derivatives of constraints to the main evolution sys- 
tem. We studied the initial-boundary value problem for 
two different formulations of Maxwell's equations, the fat 
Maxwell system and the Knapp-Walker-Baumgarte for- 
mulation. We noticed that the modification alters the 
character of the systems, transforming them into mixed 
hyperbolic-parabolic problems. Parabolicity introduces 
two remarkable effects. First, it damps the high fre- 



7 



- u 
W 



original 

modified 



10 
t 



- u 
W 10 



original 

modified 



10 

t 



io 2 



FIG. 3: At a fixed resolution of h — 1/50 and with periodic 
boundary conditions and no dependence in z we compare the 
decay of the constraint energies, Eqs. (11 21 and 1)27^ . between 
the modified and unmodified systems for the fat Maxwell 
(top) and the KWB (bottom) systems. As in Fig.|2]the initial 
data is a constraint violating pulse of compact support to all 
fields. Although this cannot be seen from the plots, at time 
t = the constraint energy of the modified and unmodified 
systems coincide. 



quency components of constraint violations (whereas in 
the KWB case the damping affects all constraints, in the 
fat Maxwell system some constraints are not damped). 
Second, it allows for Dirichlet boundary conditions for 
the constraints which naturally translate into Neumann- 
like boundary conditions for the main system. These 
boundary conditions are constraint-preserving and, if 
the evolution of the constraints can be made strongly 
parabolic, their expression is simpler than that of the 
original hyperbolic system. Furthermore, whereas other 
techniques, such as the A-system 0, require a substantial 
enlargement of the main system, this method leaves the 
number of variables unchanged. Although we were un- 
successful in obtaining an energy estimate for the main 
evolution systems, our numerical tests not only confirmed 
that the constraints arc damped, but also indicated that 



the discretized systems are stable. 

The stringent requirement on the time step k is possi- 
bly the greatest drawback of this technique. Whereas for 
hyperbolic problems k typically has to be proportional to 
the mesh spacing, h, in parabolic or mixed hyperbolic- 
parabolic problems it has to be proportional to h 2 . To 
circumvent this restriction one could use implicit or semi- 
implicit schemes. However, we have not yet experimented 
with such schemes. 

The damping terms introduced by the spatial deriva- 
tives of constraints vanish on physical (constraint sat- 
isfying) solutions. This indicates that dissipation will 
not detrimentally affect such solutions and, although con- 
straint violations can travel at an infinite speed, any con- 
straint satisfying solution will propagate information at 
a finite speed. Of course, at the numerical level, the nu- 
merical speed of propagation usually only coincides with 
the continuum one up to truncation errors and it is even 
possible for high frequency modes to travel in the wrong 
direction [T^ . 

It is important to realize that in the nonlinear or vari- 
able coefficient case, extra care is required for the success- 
ful application of this technique. The coefficients in front 
of the parabolic terms should never become negative, as 
this would immediately allow for unbounded growth and 
lead to instabilities. Artificial dissipation of the Kreiss- 
Oligcr type should still be added to the right hand side in 
order to damp high frequency components of the hyper- 
bolic sector. In addition, this technique, just like most of 
the available cures for constraint growth, cannot be ex- 
pected to cure instabilities, such as "gauge instabilities" , 
that are not constraint violating. 

The technique illustrated in this work requires a dif- 
ferent treatment of the excision boundaries, often used 
in numerical simulations of spacetimes containing black 
holes. Whereas in hyperbolic formulations no boundary 
conditions are needed at the excision surface, provided 
that this is purely outflow, parabolicity demands bound- 
ary conditions. Our expectation is that the use of Dirich- 
let boundary conditions for the constraints at the excision 
boundary will help to prevent the growth of constraints 
near this troublesome region. 

In this work we stressed how parabolicity can be ob- 
tained by adding spatial derivatives of constraints to the 
evolution equations. For systems, such as the ones stud- 
ied in this paper, that have first order constraints, this 
is the only way it can be achieved. However, for sys- 
tems containing constraints with second spatial deriva- 
tives (like the Hamiltonian constraint of General Relativ- 
ity) one may obtain parabolicity by appropriately adding 
such constraints without further differentiation. This ap- 
proach of adding constraints and/or spatial derivatives of 
constraints to the evolution equations could be pushed 
further and one could consider the possibility of adding 
higher order spatial derivatives of constraints. However, 
in general, adding second spatial derivatives of first or- 
der constraints or first spatial derivatives of second order 
constraints will not damp high frequency components. It 
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will only render the constraint propagation more disper- 
sive. On the other hand, adding third spatial derivatives 
of first order constraints or second spatial derivatives of 
second order constraints can help to dissipate the high 
frequencies. Certainly, it is preferable to avoid increas- 
ing the order of the formulation, as this would complicate 
the treatment of boundaries. 

Although the problems to which our technique was ap- 
plied do not present the full range of difficulties contained 
in Einstein's equations, based on the success of the tests 
that we presented we hope that the technique may prove 
helpful in a more complex situation. Its effectiveness in 
non-linear problems of General Relativity will be assessed 
in future works. 
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APPENDIX A: A DISCRETIZATION OF THE 
SECOND ORDER WAVE EQUATION WITH 
SOMMERFELD BOUNDARY CONDITIONS 



Consider the wave equation, 

d^{t 7 x)=d 2 Mt,x), 
on the real line with appropriate initial data, 



0(0, a?) = 
dt4>(P,x) = f {2 \x). 



(Al) 



(A2) 
(A3) 



The analysis that follows can be readily applied to the 
first order in time and second order in space formulation 



&t<l> 



T, 

die 



As an immediate consequence of integration by parts we 
have that the energy E = j{{dt4>) 2 + (d x 4>) 2 )dx is con- 
served by any solution of the wave equation. If we dis- 
cretize the second space derivative with the second or- 
der accurate centered differencing operator, obtaining the 
semi-discrete system 



-w = D + D -*» 



(A4) 



and observe that the one-sided finite difference operators 
D± satisfy the summation by parts rule 



(u,D±v)- 00 , +c 



-(D^u,v). 



where 



it follows that any solution of (|A4|) conserves the discrete 
energy 



(u,v) rta = y^UjVjh, 



1 dt ' dt 



+ (D + <t>,D + <)>) 

— oo,+oo ■ 

(A5) 



— oc,+oo 

Note that the similar expression 



E 



dt ' dt 



(A^,A)0)- 



00,+oc 



-00, +00 



which uses centered instead of one-sided differencing is 
not a conserved quantity, unless the right hand side of 
Eq. (|A4|) is replaced by D^tpj. However, this would re- 
quire a larger stencil and would not dissipate the highest 
frequency mode, 4>j = (— 1) J . 

We now introduce an artificial boundary at x = and 
take our domain to be x > 0. We are interested in finding 
approximations for the homogeneous Sommcrfcld bound- 
ary condition 



&(t,o) = Mt,o), 



(A6) 



which lead to a second order accurate scheme. The initial 
data are assumed to be compatible with the boundary 
condition, i.e., 



gn 

dx n 



0. 



x=0 



Consider the approximation 



~dt = D °*° 



(A7) 



This requires the introduction of 4>-i- However, by com- 
bining l|A7|) with the evolution equation at j = 0, 



dt 2 



we can eliminate <j>—i from the semi-discrete system. We 
obtain 



dt 2 
d 2 6n 



= D+D-c 



J = 1,2,3,. 



dt 2 

with initial data 



2 D+4> 



^■(0) = f [1) (x 3 ), 



(A8) 
(A9) 

(A10) 
(All) 



00, +00 j 



We conclude this appendix by proving that the 
semi-discrete approximation (|A8|l - (|Alip of the initial- 
boundary value problem (|A1|) - (|A3|) and (|A6|) is stable 
and second order convergent. 
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We denote by (f> e (t,x) the solution of the continuum 
problem QAl [) -l|A3 f) and (|aT>|) and fa (t) the grid fu nction 
satisfying the semi-discrete approximation (|A8|) (|A11(I . 
The norm with respect to which we wish to prove stability 
and convergence is the one associated with the scalar 
product 

ldu dv fdu dv\ 

2Mh = 2"A A" ^U' *Ji + oo ( + ' +V) °- +OC 

(A12) 

(D + <f>, D + 4>) 0i+oo 



The fact that the time derivative of 

2 



2E = ^ h , 

2 dt \dt dt 



appropriately bounded, we can estimate the norm of the 
error, 



d 2 1 dwo d 2 WQ f dw d 2 w 



dw\ 



< F 2 h 2 



l,+oo 



Integrating, we obtain the inequality 



1,+oc 



D + w,D + — 

dt ) 0, + oo 

dw Q \ 2 1 dw „ f dw 

-df) + 2^r Foh+ {-dt' F 

™\\l + \\F\\l +00 - 



l.+oo 



IS 



dt \ dt 



\\w(t)\\l < f e^F 2 {r)h 2 Ar + f e *- T ||F(r)||; >+(X> d7 
Jo •/ o 

= 0(h 4 ) 



proves that the the approximation is stable. To prove 
convergence we introduce the error grid-function, Wj (t) = 
(j>j (t) — 4> e (t, Xj ) , and observe that it satisfies 



d 2 Wj 
dt 2 


= D. 




- w j ~ 


^Fj, 3 


d 2 w 
dt 2 


-»( 




-W - 


dw \ 




= o, 




3 = 


0,1,2,... , 




= o, 




3 = 


0,1,2,... , 



where F 3 = 0{h 2 ) for j = 1,2,... and F = 0{h). Pro- 
vided that higher derivatives of the exact solution are 



showing that the scheme is second order convergent with 
respect to the norm induced by the scalar product l|A12(l . 
It is interesting to note that the simpler approximation 



D+D-<j)j , j = 1,2,3,..., (A13) 
D+fo , (A14) 



dt 2 



dt 



although stable, reduces the global accuracy of the 
scheme to first order. A similar calculation to the one 
in the proof leads to the inequality 

\\w(t)\\ 2 <0(h 2 ). 
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